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1. Introduction 

Earth-bound gamma-ray detectors make use of particle showers generated by gamma-ray in¬ 
teractions at high altitude. High-energy gamma rays (as well as hadrons and leptons) enter the 
Earth’s atmosphere and generate a cascade of secondary particles, forming an extended air shower. 
Eor gamma-rays of 0.1 — 100 TeV few air-shower particles reach ground level. However, secondary 
particles emit Cherenkov radiation resulting in illumination of a ^ 10^ m^ patch of the ground for 
a few nanoseconds. One should be aware that only a fraction of less than 10“"^ of the total shower 
energy is converted into Cherenkov photons, and quite a few of these photons get lost before hitting 
the ground. 

To reconstruct the direction and energy of the primary gamma ray and to discriminate them 
from charged cosmic rays most of the current experiments use reconstruction techniques based on 
the second moments of the distribution of pixel amplitudes in the camera [1,2]. These techniques 
are very robust and efficient, and relatively good angular resolution (< 0.1 deg at 1 TeV with 
VERITAS) can be reached with such a reconstruction technique. However, significant additional 
information can be extracted from the recorded images, resulting in improved performance. 

A more sophisticated reconstruction technique has been developed by the CAT experiment [3] 
and improved for the H.E.S.S. collaboration [4, 6]. These methods begin with the creation of a 
template library that contains the expected shower image for a given set of shower parameters. The 
template library can then be compared to the images recorded in a given event, and, by means of a 
multi-dimensional fit procedure, the best-fit shower parameters determined. 

In this paper we present the ongoing development work on an attempt to improve the perfor¬ 
mance of VERITAS, by the use of a Monte Carlo simulation based air shower model, combined 
with a ray-tracing telescope simulation. 

2. Model Generation 

This section describes the generation of the image models, i.e. the prediction of the expected 
Cherenkov light distribution in the camera focal plane for a given set of primary particle parameters. 
The production of mean shower images is divided into two steps: the generation of large Monte 
Carlo datasets and the simulation of the detector response. 

The electromagnetic air showers in the atmosphere are simulated with the program COR- 
SIKA [7]. The simulations are performed over a range of energy, zenith angle, impact distance, 
and primary interaction depth. The positions of the telescopes on the ground are on a plane taken 
perpendicular to the shower axis, this allows for the projection and the camera planes to be parallel. 
The output of the CORSIKA simulations consists of the prediction of the light distribution on the 
ground plus the arrival direction of each photon. Eor each event, the Cherenkov photons falling 
onto the mirror elements are followed individually according to their arrival times, initial direction, 
and wavelength with additional use of the atmospheric density profile, optical absorption and some 
of the detector characteristics such as its light-collecting area, photo-tube quantum efficiency. 

The images are produced for a VERITAS camera with pixels of 0.15° diameter, and the tele¬ 
scopes use a Davies-Cotton optical design. The shower images are generated for a source at the 
centre of the camera (i.e. 0° wobble-offset) and the gamma-ray longitudinal development is ori- 
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ented along the X-axis of the camera frame. Comparison with offset relative to the camera centre 
and/or off-axis events will result in a rotation and a translation in the camera frame; these transfor¬ 
mations are applied later in the fit procedure. 

The models are generated for 9 first-interaction depths, from O Xq and S-Xq (the radiation 
length in air is Xq = 36.7 g.cm“^), 90 energies, from 30 GeV to 30 TeV, and 50 impact distances, 
from 0 to 500 m. A multidimensional interpolation algorithm is used to interpolate between the 
templates, allowing production of an image template for any shower parameters within the param¬ 
eter ranges. Examples of bi-dimensional profiles of 1 TeV shower images obtained from Monte- 
Carlo simulations are shown in Fig. 1 for different values of the impact parameter. 


I Height=27km D=0m E=1TeVn | Height=27km D=100m E=1TeV | | Height=27km D=200m E=1TeV | 



Figure 1: Image template histograms for a 1 TeV primary gamma-ray at a core distance of 0 m (left), 100 m 
(centre) and 200 m (right) at first interaction depth of 27 km. The x and y axes are in degrees. The z-axis is 
given in arbitrary units. 

An additional parameter included here is the effect of the geomagnetic field on the electro¬ 
magnetic showers. We treat the geomagnetic field B at the VERITAS location, 

B = 47iiT D=10°21' / = 58°13' (2.1) 

D and 1 being the geomagnetic declination and inclination. The declination D is the angle between 
the magnetic field and the true north. The inclination 1 is the angle between the magnetic field 
vector and the horizontal plane that is tangent to the Earth’s surface at the observer’s position. 
Fig. 2 shows simulated camera images for two values of the shower direction. 

For the present, the position of the gamma-ray source within the camera field of view is ig¬ 
nored. This could be important due to the broadening of the telescopes optical point spread function 
with distance from the camera centre. 

3. Likelihood Reconstruction & Performance 

Once the full set of templates has been created, they must be compared with the observed 
images. The likelihood reconstruction performs a global fit to the telescope image data using a 
model for the expected pixel amplitude. Shower parameters are determined by a maximisation 
of an array likelihood function developed in [4]. The maximisation of this array likelihood is 
performed using a numerical non-linear optimisation technique. 

The probability density provides a model of a measured pixel signal ^ in units of photoelectrons 
(PE) given an expectation of /i, and consists of a convolution of the Poisson distribution of the 
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Height=40km D=50m E=100GeV (|)=Odeg 



I Height=40km D=50m E=100GeV ({)=180deg 



Figure 2: Image template for two values of the azimuthal angle of the shower axis direction: 0 = 0°, 
0 = 180°. The templates shown are evaluated for a gamma-ray shower with an energy of 100 GeV, an impact 
distance of 50 m, and a first interaction depth of 40 km. Note that the vertical scale (image amplitude) is 
given in arbitrary units but the same number of air showers were simulated. The geomagnetic field affects the 
development of showers. The images are diffused and distorted when the component of the field is normal 
to the shower. 


photoelectron number n, with the resolution of the photosensor. The latter is well represented by a 
Gaussian of width where Cp is the standard deviation of the detector noise and Oy is 

the width of the single PE response function [5]: 


P{s\H,C>y,(Jy)=Y, 

n 


n -p 




n\^l + 



{s-nf \ 
l{ol + na^) j 


(3.1) 


The pixel log-likelihood Lpix function is 


In Z^pix — 2 X In P (kS” I /i, Gp, Gy). 


(3.2) 


The array likelihood Larray is calculated from the sum over the pixel likelihoods in all telescopes, 

In Parray — ^ ^ “2 X In Gp, Gy) (3.3) 

tel. pix. 

This value is constructed under the assumption that if all pixels behave like independent random 
variables. The model photon reconstruction relies on the pixel-per-pixel comparison of the actual 
shower images with the ones that are predicted for a given set of parameters. With such a large 
parameter space, using an appropriate seed position for the fit is crucial to avoid getting trapped in 
a local nunimum. 

The shower reconstruction is performed in two consecutive stages. The standard moment- 
based parameters are calculated and used to define the starting point of the mininfisation. A differ¬ 
ential evolution (DE) algorithm may also be used to derive a handful of possible estimates [8, 9]. 
DE is a stochastic, population-based optimisation algorithm for solving nonlinear optimisation 
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problems. It does not require derivatives of the function and an initial population of values is 
created randomly in a region. Eventually, the members of the population converge to the point 
of highest likelihood value. The estimate that provides the best initial log-likelihood from DE or 
moment-based algorithms is used as starting point of the minimisation. 

The array likelihood must be minimised in a 6-dimensional fit over direction (xtarget^ytarget)^ 
impact parameter (xinip,yimp), primary energy, and first interaction depth Xq. Eitting is performed 
using a Levenberg-Marquardt (EM) algorithm. This is a standard technique used to solve nonlinear 
least squares problems. It is actually a combination of two minimisation methods: the gradient 
descent method and the Gauss-Newton method. The iteration continues until the convergence 
criteria or the stopping criteria are satisfied. Convergence is assumed if either the state vector does 
not change by more than a set value between iterations or the cost function does not change by 
more than a set value between iterations. We also provide a stopping criteria for when a maximum 
number of iterations of the algorithm is reached. The outputs of the minimisation procedure are the 
6 shower parameters, the correlation matrix of the fit parameters, and the uncertainties on the best 
fit parameters. 

A goodness-of-fit test addresses the question as to how well the data are described by the 
image template and the test only involves one hypothesis. To quantify the agreement between the 
measured and the expected pixel signals, we define a ;^^-like parameter 

^ _ 111 ^array ~ (ln Tarray) 

V2 X NdE 

where the average array log-likelihood for a given /i, (7p and Gy can be calculated as below 

(In Tarray) = ^ ^ / tis In Lpix X P(^|/i, C7p, Gy) 
tel. pix. 

The performance is assessed using Monte Carlo simulations. We consider events that satisfied 
the convergence criteria and no quality cuts are applied to the data. The energy bias and the energy 
resolution as function of the simulated energy and zenith angle are shown in Eig. 3. The energy 
bias is defined as the mean of the AEjE distribution. The energy resolution is defined as the 68% 
containment of the AE/E distribution. At low zenith angle, the energy bias changes little in the 
whole energy range while, at low energies and at large zenith angle, the curves show a large bias. 

The angular resolution is defined as the 68% containment radius of the reconstructed event 
positions from a point-like source. In contrast to the moment analysis, where the reduction of 
performance at large zenith angle is due to the larger fraction of nearly parallel images, the tem¬ 
plate analysis considers more information than just the image axis so it should cope better at low 
elevation. 


(3.4) 


(3.5) 
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Figure 3: Energy bias (E —Erec.)/E (top), energy resolution (middle), and angular resolution (68% con¬ 
tainment radius) shown as function of the simulated photon energy E. The energy resolution is defined as 
the 68% containment radius around the median value of the energy bias distribution. The left (resp. right) 
column shows results when the standard Hillas reconstruction technique provides (resp. does not provide) a 
starting value for the energy of the primary. 


4. Summary 

In this contribution, we present the progress of a likelihood-based reconstruction method for 
the VERITAS telescope array. The algorithm is based on an accurate pixel-by-pixel comparison of 
observed intensities with a Monte Carlo based template. This technique has been proven in the past 
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to be extremely successful, reconstructing events with more accuracy than the traditional moment 
analysis. As many fitting algorithms, the LM algorithm has some drawbacks, such as getting stuck 
in a local minimum or being dependent of the initial parameter values. To prevent these problems, 
we introduce a global minimisation algorithm (differential evolution scheme). Its performance is 
shown in Fig. 3. 

The intrinsic capabilities of the moment analysis and the template analysis (and in particular 
the hadronic rejection capabilities) can be combined together to improve the sensitivity of the anal¬ 
ysis. Since these analyses perform differently in different energy and impact parameter domains, 
more detailed studies should also allow one to select the optimal response on an event-by-event 
basis and therefore improve the quality (angular resolution,...) of the analysis. 

Additionally, the template method can also be used to reconstruct any particle type, in partic¬ 
ular iron nuclei [10]. 
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